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Engineering of synthetic magnetic flux in Bose-Einstein condensates [Lin et al., Nature 462, 628 
(2009)] has prospects for attaining the high vortex densities necessary to emulate the fractional 
quantum Hall effect. We analytically establish the hydrodynamical behaviour of a condensate in a 
uniform synthetic magnetic field, including its density and velocity profile. Importantly, we find that 
the onset of vortex nucleation observed experimentally corresponds to a dynamical instability in the 
hydrodynamical solutions and reveal other routes to instability and anticipated vortex nucleation. 
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One of the driving forces behind quantum degenerate 
gas research is the emulation of many-body condensed 
matter phenomenon pQ. For quantum degenerate Bose 
gases, achievements include the observation of Bloch os- 
cillations [2 and the Mott Insulator superfluid transi- 
tion [3 . Considerable attention has also been applied to 
the achievement of the fractional quantum Hall (FQH) 
regime in Bose-Einstein condensates (BECs) [4 j. One 
of the drawbacks to using ultra-cold atoms to emulate 
condensed matter phenomenon is the breaking of time 
reversal symmetry to emulate the effects of electromag- 
netic fields on charged particles. In the context of FQH 
physics considerable focus has been applied to the break- 
ing of time reversal symmetry through rotation. 

For BECs in rotating traps the nucleation of vortices 
has been observed by several groups [5 j. Theoretical 
studies, in the Thomas-Fermi (TF) regime, have proved 
effective in calculating the rotation frequency at which 
vortices are nucleated [BHS]- Such studies, which have 
been extended to dipolar BECs [9] , agree with numerical 
predictions from the Gross-Pittaevskii equation [8] flOl 
[T6] , However, to reach the FQH regime the number of 
vortices needs to be significantly larger than the number 
of bosons in the BEC. Typically, experiments are carried 
out in parabolic traps, defined by average in-plane trap- 
ping frequency uo_\_. When the rotation frequency, Q Z: is 
equal to uj_\_ the BEC becomes untrapped. Hence, the at- 
tainment of the FQH regime, which requires Q z — >• uj±, is 
an extremely challenging task. An alternative approach 
is to generate a synthetic vector potential, which breaks 
the time reversal symmetry of the problem, producing a 
synthetic magnetic field [T7H23] . 

The recent work of Lin et al. [23] has demonstrated 
the effectiveness of such an approach, in a 87 Rb BEC, 
through the experimental realization of a synthetic vector 
potential in the Landau gauge: A* = corresponding 
to a synthetic magnetic field B* = V x A*. The authors 
showed that vortices were nucleated at a critical synthetic 
magnetic field. Here we generalize the TF methodology 
used to successfully calculate the onset of vortex nucle- 



ation in rotating systems [7J H] to the case of synthetic 
magnetic fields in harmonically confined BECs. We an- 
alytically determine stationary solutions, including the 
aspect ratios, of a BEC and determine the critical syn- 
thetic magnetic field at which the stationary solutions 
become unstable. In the rotating case this instability 
corresponds to the onset of vortex nucleation. We find 
that this analysis predicts the synthetic magnetic field at 
which vortex nucleation occurred in the experiments of 
Lin et al [23] . 

Our starting point is a generalized form of the Gross- 
Pittaevskii equation, which provides a mean-field de- 
scription of the condensate wave function, i/j = ip(r,t): 
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where m is the atomic mass, V(r, t) is the confining po- 
tential and g = 47rh 2 a s /m defines the contact interac- 
tions, via the s-wave scattering length a s . 

Expressing x/j = y/pexp[i0] in the TF limit we can 
rewrite the Gross-Pittaevskii equation as: 
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where the generalized velocity is v 
v ft/mVO. 

To investigate stationary solutions we set dp/dt = 
dv /dt = and solve for p and v. Having found static so- 
lutions, they are not necessarily stable and so we analyze 
their dynamical stability. To do this we consider small 
perturbations {dp and SO) to the stationary solutions. 
Then, by linearizing the hydrodynamical equations, the 
dynamics of such perturbations are described by [7] E] 
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To investigate the stability of the BEC we look for eigen- 
functions and eigenvalues of the above operator: dynam- 
ical instability arises when one or more eigenvalues pos- 
sess a real positive part. The size of the real eigenvalues 
dictates the rate at which the instability grows. Imagi- 
nary eigenvalues relate to stable collective modes of the 
system, i.e., sloshing and breathing [24]. In order to find 
such eigenfunctions we consider a polynomial ansatz for 
the perturbations in the coordinates x, y and z of total 
degree n. For the examples considered below, all opera- 
tors in Eq. Q, acting on polynomials of degree n, result 
in polynomials of (at most) the same degree. Therefore, 
Eq. Q can be rewritten as a scalar matrix operator, 
acting on vectors of polynomial coefficients, for which 
finding the eigenvectors and eigenvalues is a trivial task. 

Below we consider the static solutions to Eqs. (2) and 
(3) for a BEC in a static elliptical harmonic trap un- 
der the influence of a synthetic magnetic field in the z- 
direction. The harmonic trapping potential has the form 



V(r,t) 



x — y plane the trap has a mean trap frequency uj± and 
ellipticity e. The trap strength in the axial direction is 
specified by 7 = uj z /uj_. Since the synthetic magnetic 
field is defined in terms of the synthetic vector potential 
there is more than one possible choice of A* . We consider 
the two most popular choices, the Landau and symmet- 
ric gauges. In each case we find the following form of the 
generalized velocity provides exact stationary solutions 
to Eqs. (2) and (3): v = au±(yi + xj) - <?*A*/ra. For 
the symmetric (S) and Landau (L) gauges the synthetic 
vector potential has the forms: A* = —\B z y/2 + }B z x/2 
(S) and A* = -\B z y (L), such that B* = k£*. 

Setting dp/dt = dv/dt = in Eqs. ( |2|3[ ) the stationary 
solutions for the respective gauges are: 
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and B* z = q*B*/muj ± . Eqs. ( |5|7|9| [Eqs. ( |6|8|l0| ] can 

be solved analytically to determine a in the symmetric 
[Landau] gauge as a function of e and I?*, with the solu- 
tions being independent of the axial trapping strength, 7 

[SHE]. 

In Fig. 1(a) we plot the solutions a for the symmetric 
and Landau gauges respectively, for e = 0. In each case 
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FIG. 1: (a) For a spherically symmetric trap, a as a function 
of the synthetic magnetic field strength, in the symmetric 
(solid black curve) and Landau (dashed green curve) gauges, 
(b) k XiV as a function of B*, in the symmetric and Landau 
(c) a as a function of the synthetic magnetic 



field (symmetric gauge), (d) the condensate ellipticity, (k x — 
K y)/( K x + ^y) 5 f° r the upper branch solutions in Fig. 1(c) as 
a function of the synthetic magnetic field strength [25] . In (c) 
and (d) e = 0.01 (solid black curve), e = 0.05 (dashed green 
curve) and e = 0.1 (red dotted curve). 

for B* < 2 only one solution exists with two additional 
solutions, referred to as the upper and lower branch so- 
lutions, bifurcating at B* z = 2. Despite the quantitative 
differences in the solutions for a the generalized veloc- 
ities, relevant in determining the physical properties of 
the BEC, are gauge invariant. 

The solutions for a [Fig. 1(a)] have qualitatively similar 
properties to the solutions obtained in the TF limit for a 
BEC in a rotating trap. This similarity is best observed 
through a comparison of the symmetric gauge solutions 
with solutions for a BEC in a rotatingtrap [6H5]. In a 
rotating trap a is determined by Eq. (|5|, with B*/2 re- 
placed by the dimensionless rotation frequency about the 
z-axis ( Q z = Q z /uj±), with terms proportional to (B*) 2 
in Eqs. ( |7|9[ ) set to zero. This difference emphasizes the 
key advantage a synthetic magnetic field has over phys- 
ical rotation. For physical rotation, when Cl z —> 1, uj x 
(Coy) tends to zero for the upper (lower) branch solution, 
hence the BEC is no longer confined. In the synthetic 
magnetic field case, the term proportional to (B*) 2 pre- 
vents cj Xi y tending to zero and the BEC remains confined. 
In principle, this can enable the attainment of high vor- 
tex densities and hence the emulation of FQH physics. 

We now consider how the synthetic magnetic field al- 
ters the TF condensate density profile: 
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with po = 15N/(S7rR x R y R z ). In Fig. 1(b) we plot the as- 
pect ratios of the condensate, K, 2 y = R 2 ^ y /R 2 = J 2 /& 2 ^ y 
[R 2 = 2gpo/(muj 2 )], as a function of the synthetic mag- 
netic field, for both the symmetric and Landau gauges, 
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with 7 = 1 and e = 0. The solutions for K x ^ y are gauge 
independent, with the solutions bifurcating at B z = 2. 
For low synthetic magnetic fields, B z < 2, the aspect 
ratios of the condensate in the x and y directions are 
reduced by > 25%, with Ki X K/y . For B* > 2, con- 
sistent with the solutions obtained for a we have three 
solutions for n x ^ y . Since the solutions are gauge indepen- 
dent, with no loss of generality, we describe the depen- 
dence of the aspect ratios, for B z > 2, in terms of the 
solutions for a in the symmetric gauge [Fig. 1(a), solid 
black curve]. For the a = solution the aspect ratios of 
the condensate continue to reduce with n x = n y . For the 
upper (lower) branch solutions for a k x follows the upper 
(lower) branch solutions and K y follows the lower (upper) 
branch solutions in Fig. 1(b), with K x (a) = K y (—a). 

For e 7^ the solutions for a significantly change [Fig. 
1(c)]. For the symmetric gauge, the a = solution no 
longer exists and the plot has two distinct branches. The 
upper branch (a > 0) is single valued and exists over 
all B*. For this branch the ellipticity of the BEC in the 
x — y plane is plotted as a function of the synthetic mag- 
netic field, in Fig. 1(d). We find that the BEC ellipticity 
monotonically increases from e with increasing B*. The 
lower branch (a < 0) is double valued and exists only 
when B* is greater than the backbending magnetic field: 
^(b)( e )- This lower branch backbending point shifts to 
higher B* as the ellipticity is increased [solid green in 
Fig. 2(a) [25 J. In Fig. 2(a) we have characterized the 
stability of the upper branch solutions through the evalu- 
ation of Eq. Q for polynomial perturbations up to order 
n = 5 [25]. The instability region consists of a series 
of crescents [7]. Each crescent corresponds to a single 
value of the polynomial degree n, where higher values of 
n add extra crescents, from above. At the high synthetic 
field end these crescents merge to form a main region 
of instability, characterized by large eigenvalues. At low 
synthetic magnetic fields the crescents become vanish- 
ingly small, with eigenvalues several orders of magnitude 
smaller than in the main instability region. These regions 
induce instability if they are traversed very slowly [26] . 
^z(b)( e ) an d the stability of the stationary solutions of 
the upper branch are key to understanding the response 
of the BEC to the adiabatic introduction of e or B*. 

For a fixed synthetic magnetic field as the ellipticity 
of the trap is introduced adiabatically, from zero, the 
BEC can follow two routes, depending on the value of 
B* relative to #* (b) (0). For B* < B* (b) (0) the solutions 
follow the upper branch until these solutions become dy- 
namically unstable. This route to instability is schemat- 
ically indicated by the vertical black arrow in Fig. 2(a). 
For B* z > #* (bif) (0) the BEC follows the lower branch 
from a = to negative a. However, as e is increased 
B*^(e) shifts upwards [solid green curve in Fig 2(a)]. 

When B z > B*^(e) the lower branch no longer exists. 
This route to instability is schematically indicated by the 
vertical white arrow in Fig. 2(a). For fixed ellipticity as 
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FIG. 2: (a) The logarithm of the maximum positive real eigen- 
values of Eq. Q [n = 5] for the upper branch solutions of a 
as a function of B* and e. The solid green curve plots the 
trap ellipticity (e) at which the backbending point (B^^) oc- 
curs [25] and the horizontal/ vertical arrows indicate routes to 
instability, (b) For the experimental parameters of Lin et al 
[231 127] . for the a branch which starts at a = for B* — 0, 
the maximum positive real eigenvalues of Eq. Q [n = 5] . 

the synthetic magnetic field is introduced adiabatically, 
from zero, the BEC follows the upper branch solutions 
for a, which become dynamically unstable. This route 
to instability is schematically indicated by the horizontal 
black arrow in Fig. 2(a) 

In the case of rotation [8 the analogous instabilities 
due to: (i) the dynamical instability of the upper branch 
solutions or (ii) the shifting of the backbending point 
lead to the nucleation of vortices. As such we expect 
that for: (i) the adiabatic introduction of a synthetic 
magnetic field, for a fixed trap ellipticity or (ii) the adia- 
batic introduction of trap ellipticity, for a fixed synthetic 
magnetic field, the onset of vortex nucleation can be de- 
termined by the onset of instability in the upper branch 
solutions and the evolution of the backbending point. 

To test this we consider the recent experiment of Lin 
et al. [23] where the synthetic magnetic field was gener- 
ated via an engineered, Landau gauge, vector potential 
[2T] , To engineer the vector potential an 87 Rb BEC was 
illuminated with a pair of Raman laser beams with a 
momentum difference along the x-direction, producing 
three tunable dressed states. The dressed state with the 
lowest energy corresponds to the following Hamiltonian 
along the x-direction H x = (Tik x - q* A x ) 2 / 2m* . In this 
Hamiltonian A* is the synthetic vector potential, con- 
trolled by the Zeeman shift for the atoms with synthetic 
charge q* and an effective mass for the atoms along the x- 
direction m*. To produce the desired spatially synthetic 
vector potential the Zeeman shift varies linearly in the y- 
direction, achieved by the application of a magnetic field 
gradient along y, resulting in an effective vector poten- 
tial: A x = —B*y. Including the effective mass into the 
hydrodynamical analysis: 

rn_m g*B* y \ t t^m- 

\ra* dx m* J mdy^ 

( m q*B* z y\~ c 
= — au ± y+ — )i + au±x3, (12) 

the stationary solutions are determined by 

0=(^a + B*)wl+o£jl (13) 
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with B% = g*£*/(ra*cj ± ), u% defined by Eq. rt8| and 



^ = (l + e) + ^a 2 + 4* 
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In the experiments the synthetic magnetic field is 
ramped up over 0.3s. Thus as the synthetic field is in- 
creased the solution starting at a = when B* = 
will be tracked, until it becomes dynamically unstable. 
Experimentally vortex nucleation started when the syn- 
thetic magnetic flux = 7rR x R y B*), in units of the 
flux quantum (3>o = h/q*)^ passing through the BEC 
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In Eq. ( 15 ) we have redefined the contact interactions to 



include the effective mass difference in the ^-direction, 
such that a s = a s /(2.b 1 ^ 3 l xy )^ with l xy = y/h/muj±. Us- 
ing the experimental parameters [27] the stationary so- 
lution becomes unstable when B* = 0.74 [Fig. 2(b)], 
corresponding to ^ = 11.1, in close agreement with the 
experimental observation of the onset of vortex nucle- 
ation. This strongly suggests that the observed onset of 



vortex nucleation arises from a dynamical instability as 
calculated from the stationary TF solutions [28] . 

The engineering of synthetic magnetic fields, in BECs, 
offers a new and exciting route into the emulation of FQH 
physics. The experimental advantages of such an ap- 
proach, as compared to rotating systems, is the ability 
to enter the FQH regime with the number of vortices ex- 
ceeding the number of bosons. Such an arrangement also 
provides a new route to investigate the FQH regime in 
lattice systems, without the need to rotate the lattice. 
We have shown that the nucleation of vortices, the first 
step in obtaining high density vortex systems, due to adi- 
abatic changes in the synthetic field and ellipticity of the 
trap can be calculated in the TF regime. This enables us 
to analytically determine the aspect ratios of the BEC for 
fields below the nucleation point and show that the ex- 
perimental observation of the onset of vortex nucleation 
is due to dynamical instabilities in the stationary state 
of the BEC. 
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